Wavefront sensor and method of reconstructing distorted wavefronts

ABSTRACT

A wavefront sensor includes a mask and a sensor utilized to capture a diffraction pattern generated by light incident to the mask. A reference image is captured in response to a plane wavefront incident on the mask, and another measurement image is captured in response to a distorted wavefront incident on the mask. The distorted wavefront is reconstructed based on differences between the reference image and the measurement image.

TECHNICAL FIELD

The present disclosure is related generally to wavefront sensors, and inparticular to modulated wavefront sensors.

BACKGROUND

A wavefront sensor is a device utilized to measure aberrations in anoptical wavefront. Wavefront sensors are utilized in a number ofapplications, including in fields of adaptive optics, ophthalmology,etc. For example, an adaptive optics device measures the aberrations inthe optical wavefront and then corrects the aberrations such that thelight provided to a measuring device (e.g., telescope) is collimated.

Current wavefront sensors are typically described as a Shack-Hartmanntype sensor, which is a first-order method utilized to detect wavefrontslopes, and a curvature sensing type, which is a second order method.With respect to the Shack-Hartman wavefront sensor, this type ofwavefront sensor utilizes an array of lenses that focus an incomingwavefront onto a sensor—such as a CCD array or CMOS array sensor. If thewavefront is planar (i.e., does not have any aberrations) then the spotsor points of light focused by each lenslet show no displacement, i.e.the spots locate at the center. If the wavefront is not planar, then thelocal slope of the wavefront at each particular lenslet results in aproportional displacement of the spot focused on the sensor, compared tothe spots of planar wavefront.

However, the Shack-Hartmann type wavefront sensor suffers from severaldrawbacks, including poor resolution, and poor pixel utilization. Forexample, the resolution of the Shack-Hartmann type wavefront sensor islargely dependent on the number of lenslets in the microlens array. Themaximum wavefront resolution depends on the number of lenslets in themicrolens array, and its sampling period is the inverse of the spacingbetween two neighboring lenslets, which prevents high resolutionwavefront slope detection. Therefore, it would be beneficial tointroduce a wavefront sensor that cures these deficiencies.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a block diagram of a Coded Shack-Hartmann (CSH) wavefrontsensor according to an embodiment of the present invention.

FIGS. 2 a and 2 b are block diagrams illustrating steps utilized tocalibrate and utilize the CSH wavefront sensor according to anembodiment of the present invention.

FIG. 3 illustrates the capture of a planar/reference wavefront image, adistorted/measurement wavefront image, and the difference between thereference wavefront image and the measurement wavefront image.

FIG. 4 illustrate fabrication of the mask (a binary one, in this case)utilized by the wavefront sensor according to an embodiment of thepresent invention.

FIG. 5 is a graph indicating how proximal parameter λ affectsconvergence of the algorithm according to an embodiment of the presentinvention.

FIG. 6 are graphs that illustrate the accuracy of the CSH wavefrontsensor according to an embodiment of the present invention.

DETAILED DESCRIPTION

The present invention provides a wavefront sensor utilized to measurethe phase of distorted wavefronts. In one embodiment, the wavefrontsensor—described herein as a Coded Shack-Hartmann (CSH) wavefrontsensor—utilizes a foreground mask that modulates the incoming wavefrontat diffraction scale. The resulting wavefront information is encoded inthe diffraction pattern, and is numerically decoded by a wavefrontreconstruction algorithm that relies on a simple optimization method. Inthis way, the CSH wavefront sensor provides higher pixel utilization,full sensor resolution wavefront reconstruction, and efficient parallelcomputation and real-time performance as compared with traditionalwavefront sensors such as the traditional Shack-Hartmann wavefrontsensor.

FIG. 1 is a block diagram of CSH wavefront sensor 100 according to anembodiment of the present invention. CSH wavefront sensor 100 includes amask 102, sensor 104, and processing system 106, which includesmicroprocessor 108 and memory 110. Mask 102 is located in the foregroundwith respect to sensor 104, separated by a distance z. In oneembodiment, mask 102 is a uniform random binary mask, and sensor 104 isa monochromatic sensor. In addition, the distance z is small; forexample, in one embodiment the distance z may be set equal toapproximately 1.6 millimeters (mm). However, in other embodiments,depending on the application the distance z may be selected within awide range of distances to create the desired diffraction pattern (e.g.,tenths of millimeters to centimeters).

The location of mask 102 modulates the in coming wavefront, whose localslopes are encoded as local movements of the diffraction pattern, whichis then captured by sensor 104. In the embodiment shown in FIG. 1 , auniform random binary mask is utilized to modulate the incomingwavefront. However, in other embodiments, other types of masks may beutilized such as a non-uniform random, gray scale mask, a diffuser, or awavelength-dependent mask, each capable of modulating the incominglight. For example, a non-uniform random mask may be utilized to createa specific diffraction pattern onto sensor 104. In still otherembodiments, the type and or distance z the mask 102 is placed from thesensor 104 may depend on the wavelength of the light being captured. Thecaptured image is a result of how the distorted wavefront is diffractedthrough the mask 102 and is provided to processing system 106 forprocessing. As described in more detail with respect to FIGS. 2 a and 2b , measurement of the distorted wavefront requires first capturing thediffraction pattern of a planar wavefront which is utilized as acalibration or reference image. Subsequently, the diffraction pattern ofa distorted wavefront is captured—designated as a measurement image.Given the reference image and the measurement image—the image pair, fromwhich the local slope information of the distorted wavefront isnumerically decoded.

In one embodiment, the accuracy of the CSH wavefront sensor isdependent, at least in part, on the distinctive diffraction patterngenerated among neighbors, which in turn is determined in general by themask 102 utilized, and the setup distance z as well. However, therequirement for distinctiveness over neighbors must be weighed againstpractical issues of light-efficiency and fabrication easiness. In oneembodiment, mask 102 utilizes a uniformly random pattern. However, inother embodiments the mask may be optimized to satisfy thedistinctiveness requirements while maintaining the requiredlight-efficiency and fabrication requirements. Also, in otherembodiments the mask need not be binary, for example a greyscale mask,or a colored mask that is capable for wavelength-variant amplitudemodulation, for the sake of specific applications.

In one embodiment, the distance z is dependent on specific requirement.If the distance z is too small, this results in a small diffraction ofthe mask, and thus weakens the distinctiveness between neighbors. Thisis best illustrated via a hypothetical in which the distance z isreduced to zero, in which case the captured image is always equal to themask pattern itself. In addition, if the distance z is too small, thenthe warped displacement flow z∇ϕ(r) becomes too small and isnon-sensible by the Horn-Schunck intensity preserving assumption, and ismore vulnerable to the sensor pixel noise and quantization error. Incontrast, a large z value makes the CSH wavefront sensor physicallybulky, and more diffractive is the unknown wavefront scalar field. Thisweakens our basic assumption that the unknown wavefront phase is oflow-diffraction. Further, large z value also enhances the diffraction ofmask 102, and possibly leads to a low contrast diffraction pattern andthus weakens the distinctiveness between neighbors as well, consequentlyimpair the wavefront sensing ability of the CSH.

FIGS. 2 a and 2 b are block diagrams illustrating steps utilized tocalibrate and utilize the CSH wavefront sensor according to anembodiment of the present invention. FIG. 2 a illustrates thecalibration step, in which planar wavefront 200 is provided to CSHwavefront sensor 100 to allow for the capturing of a reference image. Inone embodiment the captured reference image is stored to memory 110. Inone embodiment, the illumination source utilized to generate the planarwavefront is a collimated source. Once the reference image has beencaptured. CSH wavefront sensor 100 is ready to be utilized. An exampleof a captured reference image is illustrated in FIG. 3 a , below.

FIG. 2 b illustrates the measurement step, in which distorted wavefront202 is provided to CSH wavefront sensor 100 to allow for the capturingof a measurement image. As discussed above, mask 102 is interpreted asmodulating the incoming distorted wavefront 202. After propagatingthrough a small distance z, the wavefront is encoded and then capturedby sensor 104. Once again, in one embodiment the captured referenceimage is also stored to memory 110. Having captured the measurementimage, processing system 106 compares the diffraction pattern of thereference image with the diffraction pattern of the measurement image.As discussed in more detail below, microprocessor 108 (shown in FIG. 1 )utilized to compare the diffraction patterns may be implemented with acentral processing unit (CPU) or a graphical processing unit (GPU). TheGPU provides parallelism benefits over the CPU, and therefore mayprovide shorter execution times than the CPU. The measured differencebetween the two images results in sub-pixel shifting, which indicatesthe changes of local slopes. A benefit of this approach is that themeasurement is independent of wavelength, so white incoherentillumination is sufficient and no coherent lights are needed.

The following derivation is provided to illustrate the feasibility ofthe above-identified CSH wavefront sensor, Consider a general scalarfield u₀(r) (coordinate r=[x y]T) with wave number k=2π/λ at the maskposition:u ₀(r)=p ₀(r)ƒ₀(r),  (1)where p₀(r) is the mask amplitude function (with sharp boundaries),ƒ₀(r)=exp[jϕ(r)] is the scalar field with the distorted wavefront ϕ(r)representing the value to be measured. To simplify, the followingassumptions are relied upon:

-   -   Mask p₀(r) is of high frequency (and uniformly random binary),        whose Fourier transform P₀(ρ) (where ρ is the Fourier dual of r)        is broadband.    -   Distorted scalar field ƒ₀(r) is smooth enough such that its        spectrum is bandlimited, and decays sufficiently enough to zero        in high frequency regions. The wavefront ϕ(r) is of        low-frequency, and is first-order differentiable.

Assuming the mask 102 is placed close to sensor 104 (shown in FIG. 1 )(e.g., distance z≈3 mm), that the Huygens-Fresnel based principle isrelied upon. Using a compact form of Rayleigh-Sommerfeld diffractionformula, and expanded into the Fourier domain, the diffractive scalarfield may be expressed as follows:

$\begin{matrix}{{u_{z}(r)} = {{{\exp\left\lbrack {{jk}{z\left( {1 + \frac{\nabla^{2}}{k^{2}}} \right)}^{1/2}} \right\rbrack}{u_{0}(r)}} \approx {\int{\exp\left( {j2\pi r^{T}\rho} \right){\exp\left\lbrack {jk{z\left( {1 - {\lambda^{2}{\rho }_{2}^{2}}} \right)}^{1/2}} \right\rbrack}{P_{0}(\rho)} \times \text{ }{\exp\left( {- {jkz}} \right)}{\exp\left\lbrack {jk{z\left( {1 + \frac{\nabla^{2}}{k}} \right)}^{1/2}} \right\rbrack}{f_{0}\left( {r - {\lambda z\rho}} \right)}d\rho}}}} & (2)\end{matrix}$Equation (2) may be further simplified by ignoring the diffraction ofscalar field ƒ₀(r−λzρ). This is appropriate for the following reasons:(i) the distance z is small, and the diffraction effect is not obvious,and (ii) distorted scalar field ƒ₀ (r−λzρ) is smooth and its spectrumdecays sufficiently fast, so its Laplacian matrix exponential expansionis very small and negligible as well. Consequently, the followingapproximation can be relied upon:

$\begin{matrix}{{\exp\left\lbrack {jk{z\left( {1 + \frac{\nabla^{2}}{k}} \right)}^{1/2}} \right\rbrack} \approx {\exp({jkz})}} & (3)\end{matrix}$Subsequently, scalar field ƒ₀(r−λzρ) may be approximated by imposingTaylor expansion and linearizing ϕ(r−λzρ) around r, which yields:

$\begin{matrix}{{f_{0}\left( {r - {\lambda z\rho}} \right)} \approx {{\exp\left\lbrack {j{\phi(r)}} \right\rbrack}{{\exp\left\lbrack {{- j}2\pi\frac{z}{k}\rho^{T}{\nabla{\phi(r)}}} \right\rbrack}.}}} & (4)\end{matrix}$

With these approximations, Eq. (2) can be simplified as:

$\begin{matrix}{{u_{z}(r)} = {{\exp\left\lbrack {j{\phi(r)}} \right\rbrack}{p_{z}\left( {r - {\frac{z}{k}{\nabla{\phi(r)}}}} \right)}}} & (5)\end{matrix}$where p_(z)(r−(z/k)∇ϕ(r)) is the propagated diffractive scalar field ofa linearly distorted mask p₀(r−(z/k)∇ϕ(r)), predicted by the strictRayleigh-Sommerfeld diffraction formula. It is worth noting that asindicated by Equation (5), in dispersion-negligible materials,ϕ(r)=ko(r) thus the resultant intensity is irrelevant to wave number k,indicating the possibility of using incoherent illumination for opticalpath measurements. In following contents no distinction is made betweenoptical path and wavefront, and as a result o=1/kϕ(r)=ϕ(r).

CSH wavefront sensor is based on result provided by Eq. (5). Thereference image is captured under collimated illumination (i.e., whereinu₀(r)=p₀(r)), such that equation (5) above is reduced to:I ₀(r)=|p _(z)(r)|²  (6)which represents the diffraction pattern of the binary mask p₀(r). Themeasurement image is captured under distorted illumination ƒ₀(r), suchthat equation (5) above is reduced to:I ₀(r)=|p _(z)(r−z∇ϕ(r))|² =I ₀(r−z∇ϕ(r))  (7)

In this way, the local slope ∇ϕ(r) is encoded in the diffraction patternimage pair I₀(r) and I(r), which can be solved in a number of differentways. For example, in one embodiment decoding involves utilizing asimple least squares method linearized around r (as described below). Inanother embodiment, an iterative nonlinear least squares updating usingwarping theory. After getting the wavefront slope ∇ϕ(r), the unknownwavefront ϕ(r) can then be reconstructed by integrating over thepreviously obtained slopes.

In some embodiments, however, the distance z is unknown and thedisplacements in the diffraction images pair depend on the physicalscale of r. Therefore, under these conditions the reconstructedwavefront ϕ(r) can only be determined up to a scaled constant. In thisembodiment, to correctly retrieve the wavefront ϕ(r), there is anadditional constant that needs to be computed. One method of computingthis scaling constant is to determine empirically in a least squaressense over all the wavefront samples against the ground truth wavefront(described in more detail below). In this way, the CSH wavefront sensoris regarded as a first-order sensing methods.

To numerically decode the hidden phase information of the CSH wavefrontsensor, equation (7) is further linearized around r, as follows:I(r)−I ₀(r)≈z(∇I ₀(r)^(T)∇ϕ(r)  (8)which represents the final image formation model that contains thewanted phase information ϕ(r). In this embodiment, phase informationϕ(r) is solved for directly. A benefit of this approach is that it takesadvantage of the curl-free property of the flow, and allows for solvinga joint optimization problem simultaneously.

In one embodiment, the decoding operation formulates the wavefrontretrieval problem in discrete form. In this embodiment, let the totalsensor pixel number be M, and the unknown wavefront ϕ∈

^(N), where N>M includes unobserved boundary values. By defining imagegradient fields g_(x)∈

^(M) and g_(y)∈

^(M) as the average x and y gradients of the reference image andmeasured image, and g_(t)∈

^(M) as the temporal gradient between the two, a concatenated matrixG=[diag(g_(x))diag(g_(y))]∈

^(M×2M) is obtained, where diag(·) denotes the diagonal matrix formed bythe corresponding vector. The gradient operator is defined as ∇=[∇_(x)^(T)∇_(y) ^(T)]^(T):

^(N)→

^(2N) where ∇_(x) and ∇_(y) are discrete gradient operators along x andy directions respectively. With a binary diagonal matrix M_(o)∈

_({0,1}) ^(2M×2N) denoting the measurement positions, and define

$\begin{matrix}{{{M = {\begin{bmatrix}M_{o} & 0 \\0 & M_{o}\end{bmatrix} \in {\mathbb{R}}_{\{{0,1}\}}^{2M \times 2N}}},{{we}{have}\text{:}}}{{{\underset{\phi}{minimize}{{{{GM}{\nabla\phi}} + g_{t}}}_{2}^{2}} + {\alpha{{\nabla\phi}}_{2}^{2}}},}} & (9)\end{matrix}$where the classical Horn-Schunck intensity consistency term is imposed,with a regularization term on the phase smoothness. The parameter α>0 isa tradeoff parameter between the two terms. Exploiting proper splittingstrategy, Eq. (9) can be efficiently solved using Alternating DirectionMethod of Multipliers (ADMM).

In one embodiment, a slack variable w∈

^(2N) is introduced to allow the original object function shown in Eq.(9) to be split into two parts, represented as:

$\begin{matrix}\begin{matrix}{{{\underset{\phi,w}{minimize}f(\phi)} + {g(w)}},} \\{{{subject}{to}w} = {\nabla\phi}}\end{matrix} & (10)\end{matrix}$where the two functions ƒ(·) and g(·) are defined as:ƒL

^(N)→

ƒ(ϕ)=α∥∇ϕ∥₂ ²g:

^(2N)→

ƒ(w)=∥GM∇ϕ+g _(t)∥₂ ², with w|=∇ϕ

Using the Alternating Direction Method of Multipliers (ADMM) yieldsAlgorithm (1), where ζ∈

^(2N) is the dual variable. The updating of ϕ and w may be computedusing the ADMM as follows:

Algorithm 1. ADMM for solving Eq. (10)

-   -   1: Initialize w⁰ and ζ⁰, set λ>0 and K    -   2: for k=0, 1, . . . , K−1 do    -   3:

$\left. \phi^{k + 1}\leftarrow{{\arg\min\limits_{\phi}{f(\phi)}} + {\frac{1}{2\lambda}{{{\nabla\phi} - w^{k} + \zeta^{k}}}_{2}^{2}}} \right.$

-   -   4: w^(k+1)←prox_(λg)(∇ϕ^(k+1)λ^(k))    -   5: λ^(k+1)←λ^(k)+∇ϕ^(k+1)−w^(k−1)    -   6: return ϕ^(K)        (i) The ϕ-Update Step

In one embodiment, the ϕ-update step (wavefront estimation) involvessolving a Poisson equation, which usually requires proper boundaryconditions in conventional approaches. However, in one embodiment,because of the existence of M, the unknown boundary values (i.e. thenon-zero elements of (I−M₀ ^(T)M₀)ϕ, where I is identity matrix) areimplicitly determined by minimizing the objective. When trivial boundaryconditions are assumed, the resultant Poisson equation solution leads tonon-trivial boundary values on the observed part of wavefront ϕ, whichallows for the estimation M₀ϕ. In some embodiments, the Neumann boundarycondition suffices to yield a good estimation constrained to smallestpossible N. In these embodiments, by just assuming Neumann boundarycondition on the linear operators, denoting

_(DCT) and

_(DCT) ⁻¹ as forward and inverse Discrete Cosine Transforms (DCT)respectively, ϕ-update is given by:

$\begin{matrix}\begin{matrix}{\phi^{k + 1} = {{\arg\min\limits_{\phi}\alpha{{\nabla\phi}}_{2}^{2}\phi} + {\frac{1}{2\lambda}{{{\nabla\phi} - w^{k} + {\zeta^{k}_{2}^{2}}}}}}} \\{= {\left( {\left( {{2\lambda\alpha} + 1} \right)\nabla^{2}} \right)^{- 1}{\nabla^{T}\left( {w^{k} + \zeta^{k}} \right)}}} \\{= {\mathcal{F}_{DCT}^{- 1}\left( \frac{\mathcal{F}_{DCT}\left( {\nabla^{T}\left( {w^{k} + \zeta^{k}} \right)} \right)}{\left( {{2\lambda\alpha} + 1} \right){\mathcal{F}_{DCT}\left( \nabla^{2} \right)}} \right)}}\end{matrix} & (11)\end{matrix}$where the division is element-wise. Note that forward/inverse DCT can beefficiently implemented via forward/inverse Fast Fourier Transforms(FFT), respectively.(ii) w-Update Step

In one embodiment, the w-update step (slack variable) involves theevaluation of prox_(λg)(u), the proximal operator of g(w) with parameterλ, which is defined and computed as:

$\begin{matrix}\begin{matrix}{{{pro}{x_{\lambda g}(u)}} = {{\arg\min\limits_{w}{{{GMw} + g_{t}}}_{2}^{2}} + {\frac{1}{2\lambda}{{w - u}}_{2}^{2}}}} \\{= {\left( {1 + {2\lambda M^{T}G^{T}GM}} \right)^{- 1}\left( {u - {2\lambda M^{T}G^{T}g_{t}}} \right)}} \\{= {{{M^{T}\left( {I + {2\lambda G^{T}G}} \right)}^{- 1}\left( {{Mu} - {2\lambda G^{T}g_{t}}} \right)} + {\left( {1 - {M^{T}M}} \right)u}}}\end{matrix} & (12)\end{matrix}$wherein I+2λG^(T)G is diagonal in blocks and whose inverse is alsodiagonal in blocks and in closed-form, so is prox_(λg)(u). In thisembodiment, the computation of prox_(λg)(u) is element-wise. In oneembodiment, to suppress noise, a median filtering is imposed to thefinal gradient estimation. In one embodiment, the estimated wavefrontϕ_(estimate) is computed by solving Eq. (11), with median filtered(w^(K)+ζ^(K)). In this embodiment, M₀ϕ_(estimate) is used as the finaloutput wavefront.

One of the benefits of the described decoding process is the ability toexecute the decoding process in a highly efficient and parallelizablemanner. For example, the update step described above exploits FastFourier Transform (FFT) operations, element-wise operation, or simpleconvolutions. In addition, in one embodiment processing system 106(shown in FIG. 1 ) utilizes GPU to perform the FFT operations. In oneembodiment, the CUDA and CuFFT library is utilized to execute FFToperations, and the unknown size N is set to a power of two. Given theactual pixel number M, the unknown size N can be defined arbitrarily aslong as N>M, and hence is well suited for numerical implementations. Inone embodiment, utilization of N−M˜10 is utilized, although in otherembodiment other value of N may be utilized.

Several of the parameters described above may be modified, with eachmodification affecting the subsequent influence on the convergence speedof the decode algorithm. For example, in one embodiment the proximalparameter λ determines, at least in part, the convergence speed of theADMM algorithm shown in FIG. 5 .

FIG. 3 illustrates the capture of a planar/reference wavefront image, adistorted/measurement wavefront image, and the difference between thereference wavefront image and the measurement wavefront image. Inparticular, the images on the left illustrate a captured referenceimage, wherein the top portion illustrates the entire captured image andthe bottom portion illustrates a magnified view of the top portion. Theimages in the center illustrate a captured measurement image, whereinthe top portion illustrates the entire captured image and the bottomportion illustrates a magnified view of the top portion. The images onthe right illustrate the difference between the captured reference imageand the captured measurement image. The difference between the referenceimage and the measurement image illustrates the local shifting of thediffraction pattern as a result of the distorted wavefront. In this way,FIG. 3 illustrates the necessity of capturing a reference image based ona planar wavefront in order to prepare the wavefront sensor for capture.

FIG. 4 illustrates fabrication of a binary mask utilized by thewavefront sensor according to an embodiment of the present invention. Atstep 400, mask pattern 402 is written to a first side of photomask 404.In this embodiment, photomask 404 is a Soda Lime photomask, but in otherembodiments various other well-known types of photomasks may beutilized. In addition, in this embodiment a laser direct writer isutilized to write mask pattern 402 onto the first side of photomask 404.

At step 406, a thin layer 408 is deposited onto a side of wafer 410. Inthe embodiment shown in FIG. 4 , thin layer 408 is a chromium (Cr) layeris deposited onto a translucent substrate 410. In one embodiment,translucent substrate 410 is a fused silica wafer, but in otherembodiments other types of translucent substrates could be utilized. Inone embodiment, the translucent substrate may be completely transparent.

At step 412, a uniform photoresist layer 414 is deposited on top of thethin layer 408. In one embodiment, a spin-coating method is utilized todeposit the photoresist layer onto thin layer 408.

At step 416, wafer 410 is aligned with photomask 404 on a contactaligner and the photoresist layer 414 is exposed to UV light. Dependingon the type of photoresist layer utilized (positive or negative), eitherthe areas exposed to UV light become soluble to the photoresistdeveloper or the areas exposed to the UV light become insoluble to thephotoresist developer.

At step 418, a development solution is utilized to either remove areasexposed to UV light at step 416 or remove areas unexposed to UV light atstep 416.

At step 420, an etching process is utilized to etch portions of thinlayer 408. Having etched the exposed portions of thin layer 408,remaining photoresist is removed, leaving the final binary mask formedon fused silica wafer 410.

Accuracy of a proposed Coded Shack-Hartmann wavefront sensor is providedbased on experimental results. The experimental setup relies ongeneration of a plurality of different target wavefronts. In thisembodiment, target wavefronts are generated using a reflectiveLCoS-based Spatial Light Modulator (SLM) (e.g., PLUTO SLM by HOLOEYE)having a pixel pitch of approximately 8 μm, and the maximum phaseretardation of 2π. For collimated illumination, a faraway white pointlight source is employed, along with a polarizer placed in the opticalpath to allow for phase-only modulation of the SLM. A beamsplitter isemployed, along with a Kepler telescope configuration utilized to ensurethe SLM and wavefront sensor are in conjugate, to ensure that thedistorted wavefront measured by the CSH wavefront sensor is the oneproduced by the SLM. In this experiment, the Kepler telescope iscomposed of two lenses, with focal lengths of 100 mm and 75 mm,respectively. The sensor exposure time is set to be 5 ms in bothcalibration and measurement steps.

Using the above mentioned optical setup, four types of wavefront sampleswere generated and acquired, each with a series of different scales: (i)cubic phases, (ii) spherical phases, (iii) single-mode Zernike phases,and (iv) customized Zernike phases. Each of these wavefronts areanalyzed to determine how parameters may be tuned to measure and detectthe desired wavefront, as well as the tradeoff expected betweenexecution time and accuracy of the detected wavefront.

Tuning of Proximal Parameter λ

FIG. 5 is a graph indicating how proximal parameter λ affectsconvergence of the algorithm according to an embodiment of the presentinvention. Tuning of proximal parameter λ determines the convergencespeed of the ADMM algorithm. In one experiment, proximal parameter λ isevaluated for a fixed iteration number K=20, using the sphericalwavefront data. Results are provided in FIG. 5 , in which a plurality ofvalues of proximal parameter λ are evaluated. As illustrated in FIG. 5 ,a large value of proximal parameter λ makes the problem stiff andtherefore convergences slowly, while a smaller value of proximalparameter λ slows the convergence when iterations increase because of alack of constraint. In this experiment, a suitable choice of proximalparameter λ lies in the range of (0.002, 0.02), with a value of λ=0.01being utilized herein.

Tradeoff Between Time and Accuracy

In addition to the proximal parameter λ, there is also a tradeoffbetween time and accuracy of the decoding algorithm. For example, tosuppress noise, a median filtering on the estimated gradients is imposedat the last step. In practice, a median filtered final gradientestimation could indeed benefit for a smoother reconstructed wavefront.In addition, lower energy, in this case in the form of a large number ofiterations and more computation time, does not equivalently yield betterresults. On the other hand, the original energy function may notconverge sufficiently if given inadequate number of iterations.Consequently, it is necessary to investigate the tradeoff between timeand accuracy. In one embodiment, a value of K=10 yields a good tradeoff.

Scaling Constant

As mentioned before, a scaling constant is also computed. In oneembodiment, the scaling constant is determined by way of (i) meannormalization on the reconstructed wavefront; (ii) comparing with groundtruth and determining a suitable scaling constant for each pixelposition, (iii) getting the median as the scaling constant for eachindividual wavefront reconstruction, (iv) and obtaining the finalscaling constant as the mean of the medians in step (iii). In oneembodiment, an estimation for the scaling constant is set to 0.25.

FIG. 6 are graphs that illustrate the accuracy of the CSH wavefrontsensor according to an embodiment of the present invention. The topgraph illustrates mean absolute error of all sample wavefronts. Thebottom graph illustrates average angular error of all the samplewavefront gradients. As illustrated in the top graph, for most wavefrontsamples, the mean absolute error is less than 0.1 wavelengths, which is˜50 nm if monochromatic green light is considered. This illustrates thatthe CSH wavefront sensor is highly sensitive.

In addition to good error performance, the CSH wavefront sensor providesefficient run-time operation. In one embodiment, the run-time of thewavefront reconstruction algorithm is mainly reduced by reducing therequired number of iterations. For example, in one embodiment the firstiteration or update of the wavefront ϕ may be skipped when slackvariable w⁰ and ζ⁰ are set to zero vectors.

The present invention provides a wavefront sensor utilized to measurethe phase of distorted wavefronts. In one embodiment, the wavefrontsensor—described herein as a Coded Shack-Hartmann wavefrontsensor—utilizes in a foreground mask that is utilized to modulate theincoming wavefront at diffraction scale. The resulting wavefrontinformation is encoded by the diffraction pattern, and is numericallydecoded by a wavefront reconstruction algorithm that relies on a simpleoptimization method. In this way, the Coded Shack-Hartmann wavefrontsensor provides higher pixel utilization, full sensor resolutionwavefront reconstruction, and efficient parallel computation andreal-time performance as compared with traditional wavefront sensorssuch as the traditional Shack-Hartmann wavefront sensor.

While the invention has been described with reference to an exemplaryembodiment(s), it will be understood by those skilled in the art thatvarious changes may be made and equivalents may be substituted forelements thereof without departing from the scope of the invention. Inaddition, many modifications may be made to adapt a particular situationor material to the teachings of the invention without departing from theessential scope thereof. Therefore, it is intended that the inventionnot be limited to the particular embodiment(s) disclosed, but that theinvention will include all embodiments falling within the scope of theappended claims.

The invention claimed is:
 1. A method of measuring a phase of adistorted wavefront, the method comprising: capturing a reference imagebased on a planar wavefront directed through a foreground mask to asensor located behind the foreground mask; capturing a measurement imagebased on a distorted wavefront directed through the foreground mask tothe sensor; forming a diffraction pattern image pair with the referenceimage and the measurement image, wherein a local slope of the distortedwavefront is encoded in the diffraction image pair; and numericallysolving a joint optimization problem, including comparing the referenceimage to the measurement image to decode the local slope of thedistorted wavefront from the diffraction pattern image pair.
 2. Themethod of claim 1, wherein the reference image is captured utilizing acollimated light source.
 3. The method of claim 1, wherein the distortedwavefront is reconstructed utilizing an alternating direction method ofmultipliers (ADMM) methodology.
 4. The method of claim 1, wherein theforeground mask is selected from the group comprising a binary mask, auniform random binary mask, a gray scale mask, diffuser, and awavelength-dependent mask.
 5. The method of claim 4, wherein theforeground mask is located a distance z in front of the sensor.
 6. Themethod of claim 1, wherein the foreground mask is a uniform randompattern mask that modulates an incoming wavefront to create an encodeddiffraction pattern.
 7. The method of claim 1, wherein numericallysolving the joint optimization problem includes a least squares methodlinearized around a scalar coordinate r.
 8. A method of measuring aphase of a distorted wavefront, the method comprising: capturing areference image I₀(r) based on a planar wavefront directed through aforeground mask to a sensor located behind the foreground mask;capturing a measurement image I(r) based on a distorted wavefrontdirected through the foreground mask to the sensor; and numericallysolving a joint optimization problem using the reference image and themeasurement image to obtain a local slope ∇ϕ(r) of the distortedwavefront, including linearizing a diffractive scalar field uz(r) of thedistorted wavefront around scalar coordinate r.
 9. The method of claim8, wherein the diffractive scalar field uz(r) is represented as:${{u_{Z}(r)} = {{\exp\left\lbrack {jk{z\left( {1 + \frac{\nabla^{2}}{k^{2}}} \right)}^{\frac{1}{2}}} \right\rbrack}{u_{0}(r)}}}.$10. The method of claim 9, wherein linearizing the diffractive scalarfield uz(r) of the distorted wavefront around scalar coordinate r isrepresented as:I(r)−I ₀(r)≈z(∇I ₀(r)^(T)∇ϕ(r)).
 11. The method of claim 10, wherein thereference image I₀(r) and the measurement image I(r) form a diffractionpattern image pair [I₀(r), I(r)], and wherein the local slope ∇ϕ(r) ofthe distorted wavefront is encoded in the diffraction pattern image pair[I₀(r), I(r)].
 12. The method of claim 11, wherein the foreground maskis a uniform random pattern mask that modulates an incoming wavefront tocreate encode the diffraction pattern image pair [I₀(r), I(r)].
 13. Amethod of measuring a phase of a distorted wavefront, the methodcomprising: capturing a reference image based on a planar wavefrontdirected through a foreground mask to a sensor located behind theforeground mask; capturing a measurement image based on a distortedwavefront directed through the foreground mask to the sensor; andnumerically solving a joint optimization problem using the referenceimage and the measurement image to obtain a local slope of the distortedwavefront, wherein the foreground mask includes a uniform random patternmask.
 14. The method of claim 13, wherein the reference image and themeasurement image form a diffraction pattern image pair, wherein thelocal slope of the distorted wavefront is encoded in the diffractionimage pair by the uniform random pattern mask.
 15. The method of claim13, wherein numerically solving the joint optimization problem includesa least squares method linearized around a scalar coordinate r.